Domains

A Domain is the physical space where fields and operators will be defined, associated with a coordinate system. A domain is defined with its dimension (1, 2 or 3 dimensions domains are allowed in HySoP), some geometrical properties and some boundary conditions. Each point in the domain is identified with its space coordinates, denoted as \(x, y, z\).

At the time only box-shaped domains are available in HySoP, thanks to class Box

>>> from hysop import Box
>>> # 1D
>>> b = Box(length=[2.], origin=[-1.])
>>> # 3D
>>> dom = Box(length=[1., 3., 5.], origin=[-1., 2., 3])
>>> # Default
>>> dom = Box(dim=3)
>>> print(dom.length)
[1. 1. 1.]
>>> # coordinates of origin and 'highest' point
>>> print(dom.origin, dom.end)
[0. 0. 0.] [1. 1. 1.]

Origin denotes the coordinates of the ‘lowest’ point in the box. Default dimension is 3d, default sides length is 1. in each direction, and default origin position is 0 for each direction

Boundary conditions are set to periodic by default, but other types are available through BoxBoundaryCondition.

Parallelisation of the simulation

The simulation can be parallelised in two different manners:

• by tasks, i.e. some ‘tasks’ are defined and a group of processes will be affected to each task. Each group is autonomous and work in an asynchronous way.

• by splitting the data : the domain, and the fields defined on this domain, is splitted into sub-domains, each sub-domain is attributed to a process. Then each process works on its own subset of the data, in an asynchronous manner.

Both ways can be combined. Synchronisation is possible and processes or groups can communicate. This is detailed in the following parts.

Communicators and tasks

Please check MPI interface and tools in HySoP for communicator, tasks and other mpi-related definitions.

When created, a domain is automatically associated with a communicator and the processes of the communicator with a task. Default task is HYSOP_DEFAULT_TASK_ID for all processes, and default communicator is main_comm

>>> import mpi4py.MPI as MPI
>>> from hysop import Box
>>> from hysop.core.mpi import main_comm, main_size
>>> from hysop.constants import HYSOP_DEFAULT_TASK_ID
>>> dom = Box(dim=2)
>>> print(MPI.Comm.Compare(dom.task_comm, main_comm)==1)
True
>>> # Mind the objects are not identical
>>> print(dom.task_comm == main_comm)
False
>>> print(dom.all_tasks == set([HYSOP_DEFAULT_TASK_ID, ] * main_size))
True

A process can be affected to one and only one task and dom.comm_task is the communicator associated to the task of the current process.

Ok, let us now assume that we need to define three different tasks, in a simulation run with 8 mpi processes. The idea is to bind processes 0, 4, 5 and 6 to task ‘red’, 1, 2, 3 to task ‘green’ and 7 to task ‘yellow’, as shown in the figure below:

To do so, one need to set the optional argument ‘proc_tasks’ of domain, with a list which maps processes rank to task id. Indeed, proc_tasks[4] = 12 means that process of rank 4 (in main_comm!) is attached to task number 12.

Try the following program with 8 processes to check the result:

from hysop.core.mpi import main_comm, main_size, main_rank
from hysop.constants import HYSOP_DEFAULT_TASK_ID

RED = 4
GREEN = 1

proc_tasks = [RED, ] * main_size
proc_tasks[1:4] = [GREEN, ] * 3
dom = Box(dim=3, proc_tasks=tuple(proc_tasks))

assert dom.task_comm != main_comm
print('process of rank ', main_rank, ' with task ', dom.current_task())
print('rank/main_comm = ', main_rank, ' rank/comm_task', dom.task_comm.Get_rank())

Important remarks:

• task_comm defines a different object depending on which process you are.

• the rank of a process in main_comm may be different from its rank in comm_task

Some useful methods:

• current_task() : returns the task id of the current process

• task_on_proc() : returns task id of process number i

• is_on_task() : returns true if the current process belongs to task tid

MPI topologies

Domains may be distributed among several mpi processes, to allow parallel process of the simulation. MPI cartesian topologies are used to handle the mapping between hysop objects and mpi processes. For details about MPI implementation in hysop, check MPI interface and tools in HySoP.

The domain is splitted into N parts, N being the number of processes involved in the current task. MPI processes are arranged in topological patterns, on a one-, two- or three-dimensional grid. Each process is identified with its rank in the topology and with its coordinates in the grid. The grid dimension may be any number between 1 and the dimension of the domain. The figure below shows an example with 4 mpi processes distributed in a 2-dimensional grid. The physical domain is on the left and each process on the right owns its proper sub-domain. For example process of rank 2 (P2) handles the subdomain C, and has coordinates (1,0) in the grid.

Moreover, each process knows its neighbours in all directions.

Basic usage

From the simpler usage, HySoP will deal with topologies internally.

Intermediate usage

CartesianTopology objects are used to described this mpi grid layout (we’ll see later that it also handles the local space discretisation, i.e. some meshes local to each process). By default mpi can find the ‘best’ processes distribution, depending on how many processes are available, on the dimension of the space and on the way data are saved in memory (C or Fortran order indeed). The only requirement is to provide a spatial discretization using CartesianDiscretization. So the creation routine may be called very simply

>>> from hysop.topology.cartesian_topology import CartesianTopology
>>> from hysop.tools.parameters import CartesianDiscretization
>>> d2d = CartesianDiscretization(resolution=(32,32), default_boundaries=True)
>>> topo = CartesianTopology(dom, d2d)

Such a code executed with 4 mpi processes will create a topology similar to the one describe on the figure above. Notice the discretization argument, used to specified the space discretisation, we will come back to this in the next part.

The standard methods of a topology are

>>> # return the number of process in the topology
>>> print(topo.cart_size)
1
>>> # rank of the current process in the topology
>>> print(topo.cart_rank)
0
>>> # rank of the neighbours of the current process
>>> # in direction d
>>> d = 1
>>> print(topo.proc_neighbour_ranks[:, d])
[0 0]
>>> # id of the task owning this topology
>>> print(topo.task_id)
999

Advanced usage

In addition to CartesianDiscretization, different arguments are available to create any kind of topology, according to your needs. All of them are detailed in the example below

>>> # 3d domain
>>> dom = Box(dim=3)
>>> d3d = CartesianDiscretization(resolution=(33,33,33), default_boundaries=True)
>>> # let mpi choose the best distribution
>>> topo_0 = CartesianTopology(dom, d3d)
>>> # choose the dimension of the topology. dim = 1, 2 or 3
>>> topo_1 = CartesianTopology(dom, d3d, cart_dim=2)
>>> # choose the shape, i.e how many mpi process in each direction
>>> topo_2 = CartesianTopology(dom, d3d, shape=(3,1,2))
>>> # choose which direction must be distributed
>>> topo_3 = CartesianTopology(dom, d3d, cutdir=(True, False, True))

Remarks:

• shape, cutdir, periodicity: you must specify a value for each direction of the space, whatever is the dimension of the required topology. For example, topo_2 will be a 2d mpi topology but shape length is 3. For shape, the number of process in the mpi grid must corresponds to the number of processes in dom.comm_task. The example above must be executed with 6 processes.

• shape, cutdir and dim are exclusive : choose one and only one among them.

By default, topologies are periodic in each direction. Use is_periodic argument to customize this

>>> dom = Box(dim=2)
>>> d2d = CartesianDiscretization(resolution=(33,33), default_boundaries=True)
>>> topo_4 = CartesianTopology(dom, d2d, is_periodic=(False, False))
>>> topo_5 = CartesianTopology(dom, d2d, is_periodic=(False, True))

When executed with 6 mpi processes, topo_4 corresponds to left grid and topo_5 to right grid below

Periodic topologies may be useful to work with periodic domains or data.

Keep in mind that each process has its own memory space and must use ‘message passing’ to access to other processes data.

Use a predifined mpi communicator

In some cases, it may be necessary to use a previously built mpi communicator (from another soft, from fortran subroutines …)

>>> import mpi4py.MPI as MPI
>>> from hysop.tools.parameters import MPIParams
>>> other_comm = MPI.COMM_WORLD.Dup()
>>> mpi_params = MPIParams(comm=other_comm, task_id=dom.current_task())
>>> topo_6 = CartesianTopology(dom, d2d, mpi_params=mpi_params)

Space discretization

Since domain and data may be distributed among mpi processes, we must distinguish two ‘levels’ of discretization:

• global : to set the resolution of a mesh on the whole physical domain; done thanks to CartesianDiscretization. This does not depend on the number of mpi processes. This global resolution is just a ‘modelisation’ parameter. No data or array will be allocated using this resolution. Everything concerning memory is done at process level.

• local : all the parameters of the meshes on a given mpi process, i.e. the discretization of a sub-domain.

To summarize, you must choose a global discretisation for each topology, for which some local discretisations will be computed, depending on the number of processes and on the boundary conditions.

Global discretisation

For the moment, only regular grids are available in HySop, defined with the class CartesianDiscretization

>>> from hysop.tools.parameters import CartesianDiscretization
>>> d2d = CartesianDiscretization([65, 29])
>>> print(d2d.resolution)
[65 29]
>>> d3d = CartesianDiscretization([129, 33, 257])
>>> d3dg = CartesianDiscretization([129, 33, 257], ghosts=[1, 2, 1])

The two arguments are the resolution of the grid, i.e the number of points (not cells!) in each direction, and the number of points in the ghost layer in each direction. Resolution argument is mandatory, while ghosts is set to zero by default. This is a ‘global’ description of the required discretization. What really happens in memory for fields discretize using these objects strongly depends on the number of mpi processes used for the simulation and on the boundary condition types. This will be discussed later.

Local grids

When creating a topology, a global discretization parameter must be provided. From this discretization and the topologies parameters (mpi grid layout), local grids are computed and saved in the attribute ‘mesh’ of the topology

>>> from hysop import Box, CartesianDiscretization
>>> dom = Box(dim=2)
>>> d2d = CartesianDiscretization([65, 65], default_boundaries=True)
>>> topo = CartesianTopology(dom, d2d)

To compute the local resolution, the global resolution in each direction is divided by the number of mpi process in this direction. Remaining points are appended on the last process of the considered direction. Then, if asked, some ghost points are added on local boundaries. Definition and utility of ghost points is detailed thereafter.

To clarify notations, an example is detailed in the figure below. A periodic 1d domain of length Lx, starting at point x0 is condisered. The chosen discretization has 11 points, as shown on the top part of the figure. The points are numbered from 0 to 10, which corresponds to their global index. Creating a topology based on this discretization, with 2 points in the ghost layer and executed on 3 mpi processes results in the local meshes at the bottom of the figure. Pn stands for process of rank n in the 1d topology.

We distinguish computational points, in black on the figure, from ghost points, in red. The former are ‘real’ grid points where computation will occur. They corresponds to the points of the original (global) grid. Ghost points are some extra points, some kind of fake boundaries, local to each process. Points are identified thanks to their index, either local (numbers in green) or relative to the global mesh (numbers in black).

Important remark:

• when the domain is periodic, the last point in each direction is ignored on local meshes, to save memory. Either this point is in the ghost layer or not present at all if ghost layer equal to 0. Mind this to choose the values of the discretization parameter. A common usage consists in choosing a number of points which is a power of two, plus 1, to fit with fftw requirements.

To conclude, here is a list of the most useful attributes of mesh class

>>> from hysop import Box, CartesianDiscretization
>>> d1d = CartesianDiscretization([11], ghosts=[2], default_boundaries=True)
>>> dom = Box(dim=1)
>>> topo = CartesianTopology(dom, d1d)
>>> mesh = topo.mesh
>>> # local resolution
>>> print(mesh.local_resolution)
[15]
>>> # local resolution excluding ghost points
>>> print(mesh.compute_resolution)
[11]
>>> # coordinates of the origin of the local mesh
>>> print(mesh.local_origin)
[-0.18181818]
>>> print(mesh.local_compute_slices)
(slice(2, 13, None),)
>>> # list of local indices of computational points
>>> print(mesh.local_compute_indices)
(array([ 2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12], dtype=int32),)
>>> # list of local indices of all points
>>> print(mesh.local_indices)
(array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14],
      dtype=int32),)
>>> # coordinates of local points
>>> print(mesh.local_coords)
(array([-0.18181818, -0.09090909,  0.        ,  0.09090909,  0.18181818,
        0.27272727,  0.36363636,  0.45454545,  0.54545455,  0.63636364,
        0.72727273,  0.81818182,  0.90909091,  1.        ,  1.09090909]),)
>>> # coordinates of computational points
>>> print(mesh.local_compute_coords)
(array([0.        , 0.09090909, 0.18181818, 0.27272727, 0.36363636,
       0.45454545, 0.54545455, 0.63636364, 0.72727273, 0.81818182,
       0.90909091]),)
>>> # space step
>>> print(mesh.space_step)
[0.09090909]

local_compute_indices returns a list of python slices. For example, in our example of the figure above, on process 1, compute_index is equal to [slice(2, 5)], which means that first point has local index 2 and last point local index 4. This argument can be used to call a numpy array

>>> from hysop.tools.numpywrappers import npw
>>> # init some arrays.
>>> a = npw.zeros(topo.mesh.local_resolution)
>>> # set value only for computational points:
>>> a[mesh.local_compute_slices] = 3.

local_mesh_coords or local_compute_mesh_coords returns a tuple, with coordinates values in each direction. It can be used as argument in a function, like this

>>> from hysop import Box, CartesianDiscretization
>>> import numpy as np
>>>
>>> func = lambda x, y: np.cos(x) + np.sin(y)
>>>
>>> d2d = CartesianDiscretization([11, 11], ghosts=[2, 2], default_boundaries=True)
>>> dom = Box(dim=2)
>>> topo = CartesianTopology(dom, d2d)
>>>
>>> # apply 'func' for (x,y) of each computational point of topo
>>> # on the current process
>>> res = func(*topo.mesh.local_mesh_coords)
>>> # res is an array of shape topo.mesh.local_resolution
>>> print((topo.mesh.local_resolution == res.shape).all())
True

Ghost points

Ghost-points are a very common trick to manage the behavior of numerical methods at the boundaries of local (distributed) subdomains.

Consider the left-hand figure below, where for some numerical method with a 5 points stencil, one need to compute some value at point of index 4. To do so, values of two neighboring points in each direction are required. If the grid is distributed without ghost points, values from points of global indices 6 and 7 will be held by process 2 while computation for point 5 occurs on process 1. So, some ‘false’ boundaries are added, the ghost points, in red on the figure. That means that values of points 6 and 7 are duplicated on process 1 and 2. On the former, this points are ghost points, while on process 2 they are computational points.

For a time step, the algorithm will then be more or less:

• Update ghost points : each process receive values from its neighbours in its ghost points and send values to ghost points of its neighbours.

• Compute : each process apply numerical scheme for all its computational points.

When required, the update step is handled by the operator and so hidden from the end user. But, if needed, use CartesianDiscreteFieldGhostExchanger to explicitely update ghost points:

>>> from hysop.fields.ghost_exchangers import CartesianDiscreteFieldGhostExchanger
>>> ghe = CartesianDiscreteFieldGhostExchanger("exch", topo, (res,))
>>> ghe.exchange_ghosts()

A few remarks:

• the number of numpy arrays to be synchronized must be set at init and must corresponds to the length of the list argument of the apply method,

• arrays shapes must fit with the topology local resolution,

• ghost points implies mpi communications and impact the memory print of the numerical method. Therefore they must be used only when required.